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We compare the critical properties of the two-dimensional (2D) XY model in a transverse 
magnetic field with filling factors / = 1/3 and 2/5. To obtain a comparison with recent experiments, 
we investigate the effect of weak quenched bond disorder for / = 2/5. A finite-size scaling analysis 
of extensive Monte Carlo simulations strongly suggests that the critical exponents of the phase 
transition for / = 1/3 and for / = 2/5 with disorder, fall into the 2D Ising model universality class. 
Studying the possible domain walls in the system provides some explanations for our results. 
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The frustrated XY model provides a convenient frame- 
work to study a variety of fascinating phenomena dis- 
played by numerous physical systems. One experimental 
realization of this model is in two-dimensional arrays of 
Josephson junctions and superconducting wire networks 
A perpendicular magnetic field induces a finite 
density of circulating supercurrents, or vortices, within 
the array. The interplay of two length scales - the mean 
separation of vortices and the period of the underlying 
physical array - gives rise to a wide variety of interest- 
ing physical phenomena. Many of these effects show up 
as variations in the properties of the finite-temperature 
superconducting phase transitions at different fields. Re- 
cent and ongoing experiments have measured the critical 
exponents in superconducting arrays || , opening the op- 
portunity to do careful comparison of theory and exper- 
iment. In this Letter we examine the critical properties 
of the 2D XY model for two different values of the mag- 
netic field in the densely frustrated regime (/ ^> 0) and 
in the presence of disorder. 

The Hamiltonian of the frustrated XY model is 

n = -J2j ij cos(9 i -e j -A ij ), (i) 

where 9j is the phase on site j of a square L x L lattice 

and Aij = (2n/(po) f- A ■ til is the integral of the vector 
potential from site i to site j with <f>o being the flux quan- 
tum. The directed sum of the Aij around an elementary 
plaquette ^2 Aij = 27r/ where /, measured in the units of 
0o i is the magnetic flux penetrating each plaquette due 
to the uniformly applied field. We focus here on the cases 
/ = vl 1 with p/q = 1/3 and 2/5. 

A unit cell of the ground state fluxoid pattern for these 
/ is shown in Figure |l](a) The pattern consists of di- 
agonal stripes composed of a single line of vortices for 
f = \ and a double line of vortices for / = |. (A vortex 
is a plaquette with unit fluxoid occupation, ie. the phase 
gains 27T when going around the plaquette.) The stripes 
shown in Figure |l](a) can sit on q sub-lattices, which we 
associate with members of the Z q group. They can also 
go along either diagonal, and we associate these two op- 



tions with members of the Z2 group. A common specu- 
lation for commensurate-incommensurate transitions and 
the frustrated XY model is that the transition should be 
in the universality class of the q-state (or 2q-state) Pott's 
model. We find that this is not the case because domain 
walls between the different states vary considerably in 
both energetic and entropic factors. 

Table | lists the energy per unit length a for straight 
domain walls between the various ground states at zero 
temperature. We also numerically calculated the energy 
of domain walls that are not straight. Closed domains, 
like those seen in the simulations, of linear dimension L 
from 10 to 60 unit cells in a system of size 120x120, have 
energies that scale linearly in L to very high accuracy. 
Examining two domains as a function of their distance 
apart shows only short range forces between them which 
fit an exponentially decaying function of distance very 
well ||. This strongly indicates that we can treat the 
energy of these domains as being linear in L. The other 
type of excitation (in addition to spin waves) is vacancy- 
interstitial pairs. Such pairs have logarithmic interac- 
tions and can undergo a Kosterlitz-Thouless (KT) tran- 
sition Q. First, we focus on domain wall excitations. 

The fluxoid pattern for the two lowest energy walls at 
/ = I is shown in Figure [j]. One can see from the fig- 
ure that a shift wall can be viewed as two adjacent, or 
bound herringbone walls. For / = 5 the energy of two 
herringbone walls is less than that of a single shift wall 
and hence, the shift walls are unstable and break up into 
herringbone walls. As a result, we confine our discussion 
of the f = \ case to the herringbone walls as other walls 
should not be present at large length scales. The energy 
cost for dividing an L x L lattice into two domains sepa- 
rated by a solid-on-solid (SOS) wall stretching from one 
side of the system to the other is 

T-tsingie{z} = baL + ba^\z k - %k-\\. (2) 

k 

The height variables z k , take on integer values (6 = 3 
is the shortest length segment). The partition func- 
tion, Z = J2{ Zk } ex P( — 7~L/T) can be evaluated ei- 
ther by the transfer matrix method or recursively. 
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The interracial free energy per column M is T = 
Tln[e b ' T / T tanh(6cr/(2T))]. The zero crossing of T gives 
an estimate of the critical temperature. Plugging in the 
values for the f — \ herringbone wall gives T c = 0.19J, 
in remarkable agreement with the value T c = 0.22 J found 
in the Monte Carlo simulations described below. Being 
similar to Ising walls, herringbone walls cannot branch 
into other herringbone walls, thus the set of possible do- 
main wall configurations is similar to those in an Ising 
model. We label the fraction of the system in state (s,j) 
as m s j, where s = ±1 denotes the member of Z 2 , and 
j = 1,2,3 denotes the member of Z$. Below the tran- 
sition, one state (s,i) spans the system. On this state 
sit fluctuating domains, bounded by herringbone walls, 
of each of the states (— s, 1), (— s, 2), and (— s, 3) in equal 
numbers; so the Z3 symmetry is broken for the (s,j) 
states, but not for the (— s,j) states. As the transition 
is approached from below, the domains occupied by the 
(— s,j) states grow, with smaller domains of the (s,j) 
states within them. At the transition, the Z2 symmetry 
between the ±s states is restored and, as a result, the Z 3 
symmetry for the (s,j) states is also restored. 

The Monte Carlo simulations used a heat bath algo- 
rithm with system sizes of 20 < L < 96. We computed 
between 10 7 and 3 x 10 7 Monte Carlo steps (complete lat- 
tice updates) with the largest fraction close to the transi- 
tion temperature. Data from different temperatures was 
combined and analyzed using histogram techniques || . 

If the largest fraction of the system is in state (s,i), 
then we have three Ising order parameters, Mj — (m s ^ — 
m^ s ,j)/(m Si i + m_ S)J '), j = 1 • • ■ 3. On average, these Mj 
are the same so we just take the average as M. To cal- 
culate the m CT) j, we examine the Fourier transform of 
the vortex density p k ± at the reciprocal lattice vectors 
k± = ~(1, ±1) of the ground state vortex lattices. Start- 
ing from the definition of the Fourier transform, and us- 
ing the vortex states given above, one finds Pk±/Pg = 
m±i,i+m±i^e l2 ' R / 3 + m±i^e~ l27T ^ 3 , where p g is the mod- 
ulus in the ground state. In practice, pk± is reduced by 
small short-lived regions which don't quite match any of 
the six states. Since this effect is the same for all states, 
it cancels when calculating M . Using the real and imag- 
inary parts of pk± in addition to m ±ij, calculated 
from the direct vortex lattice as in M, we can find the 
five independent m a j. 

The transition temperature is located using Binder's 
cumulant @, U = 1 - (M 4 )/(3(Af 2 ) 2 ), shown in Fig- 
ure ||(a) . For system sizes large enough to obey finite-size 
scaling, this quantity is size independent at the critical 
point. From Fig. |we find T c = 0.2185(6) J. T c can also 
be determined from the scaling equation for the temper- 
ature at the peak of thermodynamic derivatives such as 
the susceptibility, T C (L) = T c + aL- 1 ^ . We find these 
other methods give T c in agreement with that from U. 

The exponents for the specific heat a, the order pa- 
rameter /J, the susceptibility 7, and the correlation length 
v describe the usual power-law singularities in the infi- 



nite system limit. Finite size scaling |TT[ at T c applied 
to dlnM/dK gives \jv = 1.007(25), and to the 
susceptibility \ gives 7/^ = 1.743(20), and to M gives 
f3/v = 0.142(20) (these exponents are determined from 
the slopes of the lines shown in Fig. |3|(a) and (b)). This 
is in excellent agreement with the Ising values v = 1, 
7 = |, and [3 = |. Fig. || shows the collapse of the raw 
data onto the scaling function (inset) for \- 

Two previous examinations of the f = \ case Q sug- 
gested a continuous transition but did not measure crit- 
ical exponents. Lee and Lee || claim to find separate, 
closely spaced transitions, for the breaking of Z2 and Zj,. 
One explanation for their conflicting results comes from 
the small system sizes (L < 42) used in their analysis. 
Below the transition, if the dominant state is (s,i), in 
small system sizes you often do not see all three of the 
(—s,j) states in the system at the same time. This can 
give the impression of separate transitions for small sys- 
tem sizes. This impression is enhanced by the presence 
of a shoulder in the specific heat at intermediate system 
sizes §. For the larger system sizes, we see this shoulder 
merge with the main peak and for L = 84 and L = 96 it 
is no longer discernible. 

The helicity modulus Y is the quantity most closely re- 
lated to experimental measurements ||. For / ^ 0, the 
scaling of the I-V curves found in experiments is consis- 
tent with domain wall activation processes || . The the- 
ory of Nelson and Kosterlitz for the / = case predicts 
that Y should come down in a characteristic square-root 
cusp and then jump with the universal value, IIlbTkt 1^ ■ 
However, we find an exceptionally good fit (Fig ||(c)) of 
ourdatatoF-Fo = L^ fl,v M({T -T C )L X / V ) with v = 1, 
(3 = |, and Fo = 0. Clearly, Y is strongly renormalized 
from its bare value and attempting to fit scaling relations 
for the / = case || without taking this into account 
seems questionable. We see two possible interpretations 
of our result. The first is that Y only receives contribu- 
tions from the ordered part of the lattice. So compar- 
isons with the / = case should examine Y m = Y/M. 
Y m « 0.58 at the transition implying a larger than uni- 
versal jump. Alternatively, one can say that although Y 
is brought down by the presence of fluctuations in M, 
it should still jump when it crosses the universal value, 
IksT '/tt. Extrapolating the observed behavior of Y gives 
Yl^too = a\T — T C \P. This crosses the value of the univer- 
sal jump at Tkt —T c « 10~ 6 . Although we do not see 
evidence for a jump (the best fit has Y = 0), a difference 
in transition temperatures of 10 -6 would not lead to any 
observable effects for the system sizes studied here. 

While for / = i , herringbone walls are the only stable 
walls, this is not true for / = |. For / = | it is ener- 
getically favorable for two herringbone walls to bind and 
form a shift-by-one or shift-by-three wall. Binding does, 
however, have an entropic cost. To see if these walls are 
bound we consider the following model for two SOS walls: 

H d oubie{A, z} = ^{(26cr + U||<5 Zfc: o) +ba\z k - z k -i\ 

k 
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+(2ba + u±S Zkt0 )A k + V r ({A, z})}. (3) 

Zk is the separation of the walls (zk > 0), A/, is the 
number of vertical steps the two walls take in the same 
direction in the k'th column (— oo < Afc < oo). tin and 
u± are the binding energies parallel and perpendicular to 
the wall. At this stage we take V r = 0. Summing over Afc 
leaves the partition function in the form of a transfer ma- 
trix: Z — J2{z k } Ylk 1 ■ Restricting Zfc — Zk-i to or 
±1, we can derive the eigenvalues and eigenvectors of the 
matrix T explicitly. A ground state eigenvector i/v( z ) — 
e~ MZ , where is the localization length, or typical dis- 
tance separating the lines, characterizes the bound state 
of the two lines, fj, = defines the unbinding transi- 
tion at T;,. Repeating this process numerically for the 
unrestricted case (zk — Zk-i > 0) gives Tb — 0.398J for 
the shift-by-one walls and Tb — 0.442J for the shift-by- 
three walls. The free energy for these walls crosses zero 
before they unbind. Hence, at the transition we expect 
a branching domain wall structure similar to the q > 5 
Pott's models where a first order phase transition occurs. 

In their Monte Carlo simulations, Li and Teitel p"3j ] 
observed hysteresis of the internal energy when the tem- 
perature was cycled around the transition and used this 
as an argument for a first order transition at / = |. 
The most direct indication of a first order transition is 
the presence of a free energy barrier between the ordered 
and disordered states which diverges as the system size 
increases [T^ ]. The free energy as a function of energy 
is obtained using Tl (E) = — In Pl (E) where Pl (E) is 
the probability distribution for the energy generated by 
Monte Carlo simulation of a L x L system. Figure ||(c) 
shows the growth in this barrier as the system size in- 
creases from L = 20 to 80 giving clear evidence for the 
first order nature of the transition. Also, according to 
finite size scaling, the maximum of C and y should scale 
with L d for first order phase transitions pl| . We find 
this to be the case and also obtain T c = 0.2127(2) J J5|. 

We now consider the effects of disorder on the / = | 
phase transition. Taking the couplings in the Hamilto- 
nian (|l|) as Jy = J(l + £y), the are chosen randomly 
from a Gaussian distribution with a standard deviation 
5. Due to variations of the phase differences across the 
bonds, a specific realization of random bonds may favor 
a certain sub-lattice for the ground state, creating an ef- 
fective random field. To quantify the effect, we placed 
the fluxoid configuration of the ground states down on 
10 000 separate realizations of the disorder and allowed 
the continuous degrees of freedom (the phases) to relax 
and minimize the energy. We find that the energy shifts 
from the 6 = case, and these shifts fit a Gaussian distri- 
bution with mean — 0.5<5 2 L 2 and standard deviation SL. 
The difference in energy between states which were de- 
generate in the clean system is the measure of the random 
field. This difference centers on zero and has a standard 
deviation of 0.75SL for two states related by a shift and 
0.57SL for two states with vortex rows along opposite di- 
agonals. The effect of random fields on discrete degrees 



of freedom in 2d is marginal |l5| . For d > 2 there is a 
critical randomness above which random fields cause the 
formation of domains in the ground state of size ~ (,rf- 
Aizenman and Wehr have shown that this critical ran- 
domness is zero in 2d j^]. Yet, their result does not 
preclude the possibility that £rf is so large as to be un- 
observable in a finite sized sample. Indeed, experiments 
on superconducting arrays have found apparent phase 
transitions, including scaling behavior |^] in sample sizes 
of order 1000 x 1000. In our simulations with disorder at 
5 < 0.1, all systems had a low temperature state with the 
order parameter approaching unity. We will, therefore, 
ignore the effects of random fields for 6 < 0.1 assuming 
that £_rf is larger than the sample size. 

At any coexistence point of the clean system, random 
bonds result in different regions of the system experienc- 
ing average couplings slightly above or below the criti- 
cal coupling. As a result, at any given temperature the 
system will predominantly prefer either the ordered or 
disordered state wiping out the coexistence region and 
leaving only a continuous transition |l5j ]. It has been 
conjectured Jlfj] that critical random Potts models are 
equivalent to Ising models. Kardar et al. fl7| suggested a 
possible mechanism for this effect. Their position space 
renormalization group approximation suggests that the 
probability of loop formation in the fractal interface of 
the clean system vanishes marginally at a transition dom- 
inated by random bonds. The interface may have some 
finite width due to a froth of bubbles of different phases, 
but under renormalization a linear critical interface is 
obtained and, hence, an Ising transition appears. 

The fluxoid configurations from our simulations sug- 
gest that for large enough disorder, (5 > Sf) the inter- 
face is really linear, not just in the renormalized sense. 
Sf can be estimated by placing a random potential V r 
in Eq. ||. Ignoring the terms involving Afc, one obtains 
the model for wetting in the presence of disorder, solved 
by Kardar Jil| in the continuum limit. He obtained a 
new length scale due to randomness, 1/n = 2T 3 /K5 2 
where K is the renormalized stiffness The unbinding 
transition is lowered and is now defined by the condi- 
tion fi — k = 0. As T(, decreases, it eventually hits the 
transition temperature for the first order phase transition 
observed in the clean system. At this point any branched 
domain wall structure is unstable. This is just the last 
step in a process in which the effective linear interface 
becomes narrower as disorder increases. In the vicinity 
of this "final" unbinding, the Ising-type behavior of the 
system should be readily visible at any length scale. 

We have done a Monte Carlo analysis with bond disor- 
der values of 5 = 0.05 and 0.1. To calculate the average 
value of a thermodynamic quantity, we first calculate it 
for a given realization of the disorder and then do a con- 
figurational average over 10 to 15 realizations for 8 = 0.1 
and seven realizations for S = 0.05. Figure |^(c) shows 
the free energy barrier for / = | as a function of sys- 
tem size in the for 5 = 0.05, and 0.1. For 6 = 0.05, the 
barrier first grows with system size and then levels off. 
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At 5 = 0.1 the free energy barriers are essentially zero, 
indicating a continuous transition and that the system 
sizes are large enough to apply finite size scaling. Here, 
we follow the finite-size scaling methods used in jl6| . 

Figure || shows the peak values of dhiM/dK and 
X as a function of L. The slopes of these plots give 
l/v = 1.05(12) and ^jv = 1.70(12). A similar analysis 
oidM/dK gives (l-/3)/v = 0.94(10) f|. Within errors, 
these exponents are what one would expect from an Ising 
model. Experiments at / = | || also found a continuous 
transition and measured the critical exponents v — 0.9(5) 
and the dynamic critical exponent z = 2.0(5), consistent 
with an Ising transition. 

In conclusion, we find that the nature and universality 
class of the phase transitions are quite sensitive to the 
proximity of the binding transition for the lowest energy 
domain walls. For / = 1/3 the lowest energy walls are 
never bound and the transition is Ising- like. For / = 2/5 
domain walls can lower their free energy by binding to 
each other, resulting in a first order phase transition. Dis- 
order weakens this binding and changes the transition to 
be continuous and Ising-like. Our results are consistent 
with the continuous phase transition and critical expo- 
nents observed experimentally for / = 2/5 S. 

We thank M. Aizenman, P. Chandra, J.M. Kosterlitz, 
X.S. Ling, and D. Huse and for useful discussions. 
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FIG. 1. Fluxoid pattern for (a) unit cells of / = | and 
/ = §, and domain walls for f ~ \ (b) herringbone wall, 
(c) shift-by-one wall, and (d) shift-by-one wall branching into 
two herringbone walls (a vortex is shown as a dark square) . 
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FIG. 2. / = 1/3 (a) Binder's cumulant U vs T for L = 36 

to L — 84 (smaller L shown as dotted lines), and (b) x vs T 

for L = 36 to L = 84 and scaling collapse of this data (inset) 

where x = (T — T c )L 1/v , y = xL~ l/v , v = 1, and 7 = \. 
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FIG. 3. Finite size scaling plots for / = § (triangles) and 
/ = |,<5 = 0.1 (pentagons): (a) logarithmic derivative of M 
vs L, (b) x vs L. (c) Scaling collapse of Y (raw data in inset) 
where x = (T — T c )L 1/v , y = Yl/I" , v = 1, and /3 = |. 
(d) Free energy barrier vs system size for / = | and 5 = 
(squares), S — 0.05 (circles) and <5 = 0.10 (triangles). 



domain wall type 



herringbone 
shift-by-one 
shift-by-two 
shift-by-three 
shift-by-four 



energy per unit length 
/ = l/3 / = 2/5 



0.056737424 J 
0.114199976 J 
0.166666666 J 



0.086117262 J 
0.158899286 J 
0.166122315 J 
0.147648594 J 
0.198688789 J 



TABLE I. Domain wall energies. 
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